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Abstract 

Background: The calculation of diffusion-controlled ligand binding rates is 
important for understanding enzyme mechanisms as well as designing enzyme 
inhibitors. We demonstrate the accuracy and effectiveness of a Lagrangian 
particle-based method, smoothed particle hydrodynamics (SPH), to study 
diffusion in biomolecular systems by numerically solving the time-dependent 
Smoluchowski equation for continuum diffusion. 

Results: The numerical method is first verified in simple systems and then 
applied to the calculation of ligand binding to an acetylcholinesterase monomer. 
Unlike previous studies, a reactive Robin boundary condition (BC), rather than 
the absolute absorbing (Dirichlet) boundary condition, is considered on the 
reactive boundaries. This new boundary condition treatment allows for the 
analysis of enzymes with “imperfect" reaction rates. Rates for inhibitor binding to 
mAChE are calculated at various ionic strengths and compared with experiment 
and other numerical methods. We find that imposition of the Robin BC improves 
agreement between calculated and experimental reaction rates. 

Conclusions: Although this initial application focuses on a single monomer 
system, our new method provides a framework to explore broader applications of 
SPH in larger-scale biomolecular complexes by taking advantage of its 
Lagrangian particle-based nature. 

Keywords: diffusion; Smoluchowski equation; Smothed particle hydrodynamics; 
Protein-ligand interactions; Binding rates; Acetylcholinesterase 


1 Background 

In the “perfect” enzyme [1] acetylcholinesterase (AChE), the rate-limiting step for 
catalysis is diffusional encounter [2, 3]. Specifically, the active site lies at the bot¬ 
tom of a 20 A-deep gorge, and the diffusion of substrate into it is accelerated by 
electrostatic steering [4, 5], Its diffusion-limited behavior, complex geometry, and 
strong electrostatic influence has made AChE a useful target for both experimental 
and computational studies of biomolecular diffusion [6, 7, 4, 5, 8, 9, 10, 11]. 

Two major classes of methods have been used to estimate diffusion rates 
in biomolecular systems. Mesoscopic coarse-grained methods like Monte Carlo 
[12, 13, 14], Brownian Dynamics (BD) [15, 8, 9], and Langevin Dynamics [16, 17] 
simulations trace the trajectories of individual coarse-grained particles driven by 
Brownian motion. Such simulations typically consider dilute ligand concentrations 
so that electrostatic protein-ligand interactions can be modeled by the Poisson- 






Pan et al. 


Page 2 of 19 


Boltzmann equation [18, 19] with a few notable exceptions [20]. Alternatively, con¬ 
tinuum models can be used to treat the diffusion of ligand concentration in space 
around a biomolecule by the Smoluchowski equation [21, 22, 23, 6, 24, 25]. In par¬ 
ticular, an adaptive finite element approach [26] has been used to numerically solve 
the Smoluchowski equation, and it shows higher accuracy in predicting experimental 
data about the ligand binding rates than the coarse-grained BD modeling [6]. For 
dilute ligand concentrations, electrostatic interactions can also be modeled with the 
Poisson-Boltzmann equation like the mesoscale approach [6, 7]. However, for more 
concentrated ligand solutions, continuum models can also model the electrostatic 
potential near the biomolecular surface using a regularized Poisson-Nernst-Planck 
formulation [24, 25], allowing screening of the ligand-receptor interactions by its 
time-dependent distribution around the protein. 

Here, we follow the continuum approach but solve the Smoluchowski equation 
using a new smoothed particle hydrodynamics (SPH) method [27, 28]. Unlike Eule- 
rian grid-based methods such as FEM, SPH is a Lagrangian particle-based method. 
SPH has been used with good accuracy for numerically solving partial differential 
equations (PDEs) describing momentum, mass and energy conservation laws [27]. 
In SPH, the domain is discretized into a set of “particles” that serve as interpolation 
points to numerically solve the governing PDEs. The SPH discretization of PDEs 
is based on a meshless interpolation scheme, which allows the PDEs to be written 
in the form of a system of ordinary differential equations (ODEs). Unlike grid- 
based FEM methods, SPH has a straightforward discretization without the need 
for time-consuming FEM mesh construction around complicated geometries such as 
biomolecules. Due to its Lagrangian nature, SPH has many advantages for modeling 
physical phenomena involving moving boundaries, large deformation of materials, 
multiphases, and advection-dominated diffusive transport [29, 30, 28]. In addition, 
the similarity of SPH to molecular dynamics and mesoscopic coarse-grained par¬ 
ticle methods (e.g., dissipative particle dynamics, BD, and Langevin dynamics), 
allows coupling of simulations across scales to build a multiscale modeling frame¬ 
work. This is our primary goal with the current work: to enable the multiscale and 
multiphysics description of biomolecular dynamics and ligand recognition. To the 
best of our knowledge, SPH has not been widely used in modeling biomolecular 
systems. Thus, in the present work, we aim to take the first step to introduce SPH 
into this Held through the development of a SPH model for biomolecular diffusion 
with AchE as a test case. 

In the SPH model, the Smoluchowski equation is numerically solved and the ligand 
binding rates are calculated from flux across the reactive boundary as in previous 
studies using FEM [21, 22, 23, 6, 24, 25] . However, in the previous FEM studies, 
active sites were modeled using the absolute absorbing (Dirichlet) boundary condi¬ 
tion (BC). This BC has a simple description on the reactive boundaries but assumes 
infinitely fast chemical reactions between the enzyme and the ligand; i.e., a “per¬ 
fect enzyme”. In our model, we take into account imperfect and non-instantaneous 
reactivity and thus solve the equation using reactive (Robin) boundary condition. 

To solve the Smoluchowski equation subject to Robin BC using SPH, we use a 
continuum surface reaction method [31] which we have recently adapted to solve the 
Navier-Stokes equations subject to slip (Robin) boundary conditions [32]. In this 
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formulation, the Robin BC is replaced by a reflective Neumann BC and a source 
term is added into the governing equation. The derivation of the method is based on 
the approximation of the sharp boundary with a diffuse interface of finite thickness 
by means of a color function. This method is general for any arbitrary complex 
geometries and thus appropriate for modeling Robin BC in biomolecular systems 
with complex structures. 

2 Results and discussion 

2.1 Spherical test systems 

Before the numerical method was applied to a biomolecular system with complicated 
geometry, we verified it on simple spherical test cases. Specifically, we considered 
a diffusing sphere with a radius R\. The entire domain was confined by the outer 
boundary Tb determined as a spherical surface with the radius of R 2 = 125 A. 

For the first test case, we let Ri = 0 and assume no external potential, for which 
the time-dependent analytical solution of the Smoluchowski equation can be easily 
derived. Figure 1 compares the SPH numerical solutions with the analytical solution 
at different times. SPH solutions are compared at different resolutions and their 
corresponding L 2 errors are calculated relative to the analytical solution. Figure 1 
shows that, even at the coarsest resolution (Ax = 8 A), the SPH solution agrees 
well with the analytical solution with about 3% relative error. This relative error is 
further reduced to 1% by increasing the resolution to Ax = 2 A. 

Next, the spherical system is assumed to have a Coulombic form of the PMF, i.e., 
W(r) = q/(3r with +1 e charge. We set R\ = 50 A and impose either a Dirichlet BC 
as specified in Eq. 6 or Robin BC as in Eq. 5. In these two tests, the corresponding 
SPH solutions of concentration at steady-state are compared with the analytical 
solution. The converged SPH solutions are shown for the Dirichlet BC (Fig. 2) and 
Robin BC (Fig. 3) imposed on the inner spherical boundary (r = R\). The reactive 
coefficient for the Robin BC is a = 1 x 10 3 . In both tests, the SPH solutions show 
very good agreement with the analytical solution even at the resolution of Ax = 8 
A, which can be further improved with increasing resolution to A.t = 2 A. Moreover, 
at Aa; = 2A, the calculated reaction rate is 2.83 x 10 12 M _1 min _1 for the Dirichlet 
BC, and is 8.24 x 10 11 M _1 min _1 for the Robin BC, both with L 2 errors less than 
3% relative to the analytically evaluated ones. 

2.2 Application to acetylcholinesterase-ligand binding rates 

We applied the SPH method to study the ligand binding kinetics of a simple spher¬ 
ical cationic ligand to acetylcholinesterase (AChE) under various ionic strength 
conditions. Specifically, we performed the time-dependent calculations at ionic 
strengths of 0.0, 0.05, 0.10, 0.15, 0.20, 0.50 and 0.67 M until the diffusion reaches 
the steady-state. To achieve the highest accuracy with affordable computational 
cost, a resolution of Aa; = 2 A was used in all the following calculations. 

In previous studies by Song et al. [6], a simple but realistic set of boundaries was 
used inspired by Tara et al. [9], encompassing the active site as well as the gorge and 
the peripheral anionic site (PAS) of AChE. We constructed these spherical active 
boundaries (T a ) at varying distances from the active site along an axis defined by 
the carbonyl carbon of S203 at the origin and the gorge. Spheres 1-6 were placed at 
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16.6, 13.6, 10.6, 7.6, 4.6, and 1.6 A along the this axis, respectively. The outermost 
spheres 1 and 2 were assigned radii of 12 and 9 A, respectively, while all others 
were given radii of 6 A. Figure 4A shows the discretized domain with i ?2 = 128 A. 
Figure 4B and 4C depict the constructed reactive boundaries 1 and 4. 

In most prior studies [6, 7, 23], an absolute absorbing (Dirichlet) boundary con¬ 
dition (Eq. 6) was assumed. However, in the present work, we demonstrated im¬ 
proved performance with the reactive (Robin) boundary condition (Eq. 5) imposed 
on the reactive boundaries. Figure 5 shows the steady-state spatial distribution 
of ligand throughout the simulation domain at different ionic strengths. At zero 
ionic strength, there are three large ligand-attracting regions, two on either side 
of the active site and one on the opposite side of the protein. There is also one 
ligand-depleted region at the top and another one near the opening of the gorge. 
At non-zero ionic strengths, electrostatic screening reduces the size of the ligand- 
enriched and ligand-depleted regions. However, a large region around the active site 
remains ligand-depleted at up to 0.50 M ionic strength. Figure 6 illustrates the tem¬ 
poral evolution of the concentration distribution as ligand moves inward in the bulk 
region from the outer boundary (T;,). The distribution has clearly reached steady 
state by 190 ns. 

We calculated the reaction rates from these solutions according to Eq. 21. In 
Figure 7, the left panel shows the time evolution of k on (t) on reactive boundary 
1 at different ionic strengths. For this boundary, k on (t) converges within 150 ns 
for all ionic strengths. The right panel shows k on (t) on reactive boundaries 1-4, 
respectively, at 0.15 M ionic strength. 

We have quantitatively compared the reaction rates calculated by SPH with ex¬ 
periment [4] and previous computational studies by FEM [8, 6]. Radic et al. [4] fit 
their experimentally measured reaction rates as a function of ionic strength using 
the Debye-Huckel limiting law 

k — (k° — k U 11 0~ 11s \ z eZi |vT lH 

"■on — V^on ' t on7 lu ' "om V 1 ! 

where I is the ionic strength, k® n is the effective reaction rate at zero ionic strength 
rate, k^ n is the effective limiting reaction rate at infinite ionic strength and set to the 
value of k on calculated at 0.67 M ionic strength, Ze is the effective enzyme charge, 
and zi is the effective inhibitor charge with a fixed value of +1 e. For the Robin 
BC SPH calculations, the reaction coefficient a was varied, as shown in Figure 8, to 
identify the value of 8.0 x 10 3 which optimized agreement between computational 
and experimental results. 

Figure 9 and Table 1 compare the reaction rates from SPH, FEM [8, 6], BD, 
and experimental data by Radic et al. [4]. As noted by Song et al. (2004) [6], 
BD simulations systematically overestimate the experimental k on , while the FEM 
produces good agreement with experimental k on at RMSD = 0.37 M _1 min _1 . With 
the Dirichlet BC, SPH predicts k on with RMSD of 0.57 M _1 min _1 , intermediate 
between FEM and BD results. However, with the Robin BC, SPH predicts k on with 
RMSD of 0.33 M _1 min _1 , better than the FEM and BD results. 

We also assessed the accuracy of SPH method for describing the ligand-binding 
kinetics of a mAchE surface mutant. We tested the surface hexa-mutant (E84Q, 
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E91Q, D280V, D283N, E292Q, and D372N) from Radio et al. [4], which reduces the 
reaction rate by about a factor of 4 across the 0 to 0.67 M ionic strengths. For the 
mutants, which are nearly isosteric with the wild-type protein, we used the same 
SPH model as the wild type, but recalculated the electrostatic potentials for the 
mutant charge distribution. As presented in Table 1, the Robin BC SPH model has 
qualitative accuracy: predicting of 2.01 compared to 1.80 from Radic et al. [4] 
with a k on RMSD of 0.70 over the entire ionic strength range studied. 

3 Conclusions 

The Robin BC offers a new way to incorporate reactive surfaces into continuum 
diffusion models for rate calculations. This Robin-based model incorporates a new 
parameter a, which has units of A//zs and can be related to the probability of 
reaction within distance Aa; of the boundary and time interval At by P = 1 — 
exp(— [33]. Thus, a = 0 corresponds to zero reactivity (reflective Neumann 
boundary conditions) while a = oo corresponds to absolute reactivity (absorbing 
Dirichlet boundary conditions). 

There are two possible origins for the differences between the current SPH model 
results and past FEM calculations using the Dirichlet BC. First, the current SPH 
work uses a more recent AChE structure (4B82) while the previous FEM calcula¬ 
tions used an older structure (1MAH). Second, our SPH model uses a fixed resolu¬ 
tion uniformly on both solution domain and boundaries, while the FEM adaptively 
meshes the reactive boundary with higher resolution. 

This work has provided an initial demonstration that the Lagrangian (particle- 
based) SPH method out-performs the Eulerian (grid-based) finite element method 
[6] in accurately predicting ligand binding rates in AChE. This result is important 
because while both methods can be used to study molecules of the size of AChE, 
SPH is more scalable to larger systems such as the synapse geometry where AChE 
operates. Additionally, due to its Lagrangian nature, SPH can easily incorporate 
other physical phenomena such as fluid flow or protein flexibility. 

We have demonstrated that superior performance can be achieved using a proba¬ 
bilistic reactive (Robin) boundary condition rather than a simple Dirichlet BC. In 
fact, the Robin BC is likely more biologically relevant than the Dirichlet BC. While 
the AChE enzyme is considered nearly “perfect” with a diffusion-limited reaction 
rate, there is experimental evidence that a very small fraction of substrates entering 
the active site gorge do not react. Specifically, recent kinetic experiments suggest 
that through unknown mechanisms, the PAS limits the rate of progression of non¬ 
substrates of any size to the catalytic site [34]. In addition, molecular dynamics 
simulations suggest that the peripheral anionic site (PAS) provides a selective gat¬ 
ing function, for example by fluctuations in the gorge width that are likely to let 
acetylcholine but not let larger molecule pass through [35, 36]. 


4 Methods 

4.1 Governing equations and boundary condition 

The time-dependent Smoluchowski equation can be written as: 


dp(x, t) 


= V • J (x, f), 


x £ fi, 


dt 


( 2 ) 



Pan et al. 


Page 6 of 19 


where p(x, t) is the concentration distribution of the reactants, and the concentra- 
tion flux J(x, f) is defined as: 

J (x, t) = -D(x) [Vp(x, t) + /3p(x, t)Vf(x)], (3) 

where D(x) is the diffusion coefficient; for simplicity, this is assumed to be constant. 
/? = 1/kfjT is the inverse Boltzmann energy with the Boltzmann constant kg and 
kinetic temperature T. W(x) is the potential mean force (PMF) for the diffusing 
particle due to solvent-mediated interactions with the target molecule. The equa¬ 
tion is solved in a three-dimensional domain D, subject to the following boundary 
conditions. First, 

p(x, t) = pbuik for x e F b , (4) 

specifying a Dirichlet BC on the outer boundary F 5 where the concentration is equal 
to a bulk concentration pbuik- The outer boundary is often a spherical surface with 
a radius chosen to ensure that the ligand-protein potential is spherically symmetric 
and/or can be approximated analytically [ 6 ]. For the current study with AChE, this 
outer boundary has radius R 2 ~ 128A as determined following a procedure similar 
to Song et al and Chen et al [ 6 , 23]. Also following Song et al and Chen et al, p is 
normalized such that pbuik = 1 - 

The active site boundary T a was modeling using either reactive Robin or absolute 
absorbing Dirichlet BCs: 

n(x) • J(x, t) = ap(x, t) for x £ T a , (5) 

or 

p(x) = 0 for x e r a , ( 6 ) 

respectively. The coefficient a is chosen to model an intrinsic reaction rate for the 
active site. Finally, a reflective Neumann BC is defined on the non-reactive boundary 
of molecule 

n(x) • J(x, t) =0 for x e T m . (7) 

Figure 10 shows the simulation domain along with all boundaries. 

Given a solution to Eq. 2, the reaction rate is calculated from the integral of the 
flux across the reactive surface [37]: 

fc on =Pbuik JI n ( x 's) • J(x' s ,t)dx' s . ( 8 ) 

r„ 

In order to solve the Smoluchowski equation (Eq. 2) subject to the reactive Robin 
BC (Eq. 5), the simulation domain is extended to include a sub-domain fl a that is 
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separated from Q by r o , and we then reformulate Eq. 2 as: 
dp r (x,t) 


dt 


=V • (D(x)[Vp r (x,t) +/3p r (x,t)VW(x)]) 

—ap r (x 1 t) JJj (n(x) + n(x')] • V x w(x — x', h r )dx!, xe!l, (9) 

f2 a 


subject to the reflective Neumann BC: 

n(x) • J r (x, t) = 0 for x £ T a . (10) 

The derivation of Eq. 9 is detailed in Appendix A, which demonstrates 


lim p r (x,t) =p(x,t). (11) 

h r —>-0 

In Eq. 9, the normalized kernel function, w(x), is a positive bell-shaped function 
with at least first continuous derivative and compact support tih r , ic(|r| > nh r ) = 0. 
The value of n depends on the specific functional form of w(x), which is specified 
in Section 4.2. In particular, w(x) satisfies the following conditions: 



fiun a 


and 


( 12 ) 


lim w(x — x', h r ) = 5(x — x'). (13) 

h r —> 0 

The normal unit vector n in Eq. 9 can be found in terms of a smoothed color 
function (j) as defined in Appendix A: 


n(x) 


V0(x) 

|V<Mx)f 


x £ U fl a . 


(14) 


4.2 SPH discretization of equations and boundary conditions 

In this section, we present SPH discretization of the Smoluchowski equation, using 
Eq. 2 if the Dirichlet BC is used and Eq. 9 if the Robin BC is assumed. To sim¬ 
plify notation, we omit superscript r for the variables in Eq. 9 in the subsequent 
derivations. 

The domain and the boundaries T a and Tf, (extended as domains Q a and Qb 
respectively), are discretized with a set of N points with positions denoted by a 
vector [i = 1, The points (which are commonly referred to as particles in 

SPH) are used to discretize and solve the governing equation. Initially, the particles 
are distributed uniformly (e.g., placed on a regular cubic lattice) with di as the pre¬ 
scribed number density at r^. The discretization is based on a meshless interpolation 
scheme: 


a.* E 

j 


—r- w(rij, h), 

a j 


( 15 ) 
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where Ai = A(ji) is a function defined at particle i , r^ and w(rij,h) is 

the weighting kernel function. The interpolation scheme assumes a summation over 
all neighboring SPH particles but, due to the compact support of w, only particles 
within distance nh from have a non-zero contribution to the summation. Spatial 
derivatives of A can be calculated as 

VjAj « ^2 h). (16) 

3 0 

In the present work, we use a cubic spline kernel as the weighting function 
f 1 - |g 2 + |g 3 0 < q < 1 

w{r,h) = ^3 | \{2-qf l<q<2 (17) 

[o q> 2, 

where q = |r| fh. With this form of weighting function, only particles within 2 h 
distance from particle i contribute to the summations in the SPH equations. We 
have chosen h = 1.3Ax where Ax is the size of the cubic lattice. 

The SPH approximation of functions and their spatial derivatives allows the 
Smoluchowski equation subject to the Dirichlet BC (Eq. 2) to be written as a 
ODE governing the evolution of concentration on particle i as: 


dpi 

dt 


E 

.jenuObuna 


A 


Dj , ,1 dw(r. 

— {Pi-Pj ) 


5 


h) 


■<*£ 
je n 


r ij 

DiPi 


drij 

D jPj 


dj 


Ti n UjT i n 


The derivations of the first and second terms on the right-hand side of Eq. 18 can 
be found in Monaghan et al. [27] where i\j is the magnitude of the vector r i3 . If the 
reactive Robin BC on reactive boundary is imposed, Eq. 9 is then solved instead 
and its corresponding SPH discretization form is: 


dpi ^ Di + Dj, , 1 dw(rij, h) 


dt d 7 - 

.jenuOt J 


: {Pi~Pj) 


rij dr.ij 


■*E 

jen 


D iPi + D jPj_^ Wi _ w _y l dw(nj,h) 


dj 


r^ drij 


- OLPi ^2 


fee £ 2 a 


+ n k 
dk 


r ik dw(r ik ,h) 
r ik dv ik 


The last term in Eq. 19 is obtained by discretizing the integral in Eq. 9 as a Riemann 


:p(x, t) [n(x) + n(x')] ’ V x w(x — x 7 , h r )dx' 


= ap{x,t) ^2 AT4[n(x) + n fe ] ■ V x w(x - r fc ,/i r ), 


(19) 


sum: 
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where A14 = ^- is the volume of particle k and Y2ken th e summation over the 
reactive boundary particles within 2 h r distance from particle i. The SPH expression 
for calculating the normal unit vector is obtained as: 


E Jr {<Pj - (/>i)Viw(r i:j , h r ) 

jGClUQa 3 

jenufi 0 3 


( 20 ) 


In the simulations presented below, we set h r = h but it could be set differently in 
future applications. 

Note that the reflective Neumann BC (Eq. 7 or 10) can be simply enforced in 
SPH by excluding the contribution from the boundary particles in the summation. 
The Dirichlet BC (Eq. 4 or 6) is enforced by assigning the fixed boundary value of 
concentration on the boundary particles. If the Robin BC is imposed, the reaction 
rate k on (t) can be calculated by Eq. 41 in Appendix B and its corresponding SPH 
discretization form is: 


k 


on 


E a Pi 

iefi 


E 

_fcer Q 


n, + n fc 
dk 


Y ik dw{r ik ,h) 
Tjj- di'ik 


( 21 ) 


Otherwise, when the Dirichlet BC is enforced on the reactive boundary, the dis¬ 
cretization of Eq. 40 in Appendix B is: 


k 


on 


E( n > • J *) 

ie fi 


E 

-fcer Q 


Yij + n fc 
dk 


Yik dw(r ik ,h ) 
Yik dr.ik 


( 22 ) 


where 


J,: = Di 


jenun b un a 


Pj - p.j Yjj dw{rjj,h) 
ru dr a 


+PD, PI y Wl-Wi luMni.h) 

< ^ r\ ■ r 


je n 


dri 


(23) 


4.3 Calculation of potentials of mean force 

We calculated the potential of mean force W(x) using the recently published 2.1 A 
resolution structure of mouse AChE [38]. To prepare this structure for the calcula¬ 
tion, we assigned titration states of ionizable residues using PROPKA [39] at pH 7, 
and we used PDB2PQR vl.8 [40, 41] to assign atomic radii and charges. APBS vl.4 
was used to perform a nonlinear Poisson-Boltzmann multi-grid calculation of the 
electrostatic potential over the entire simulation domain [42], The small and large 
domains were set to 600 A and 400 A, respectively, with a fine grid spacing of 0.600 
A. For APBS calculations, we used the single Debye-Hiickel boundary condition, 
a smoothed molecular surface, and protein and solvent dielectrics of 2 and 78.54, 
respectively. Atomic charges were mapped onto the grids using cubic B-spline dis¬ 
cretization. The calculated potential was mapped onto the SPH discretization points 
of protein and solvent via trilinear interpolation. 
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Appendix A Continuum surface reaction method 

In this appendix, we present a detailed derivation of the continuum surface reaction method for solving the 
Smoluchowski equation subject to Robin BC. We start from a two-sided problem; i.e., the concentration field p(x, t ) 
is extended into the sub-domain £l a that is separated from f2 by T a such that Eq. 2 can be approximated as 

dp (x, t ) _ ^ £) _j_ ^p r ( Xj t)VVy(x)]) — Pft (x, t) for x £ U Q a (24) 


subject to 


n(x s ) • [J r (x s ,t)| XsGr F 


j r (x s , t)l* a€r s] — o f° r x s c r a , 


(25) 


where and rf are the two sides of T a , respectively. The boundary condition Eq. 25 emphasizes that the 
extended concentration field is continuous across T a . Comparison of the weak formulations of Eq. 2 subject to 
Eq. 5 and Eq. 24 subject to Eq. 25 yields the relationship 


/// Pfj(x,t)dx = 

f2Uf2 a 



(26) 


This weak formulation is obtained by integrating the momentum equations over their respected domains and then 
applying Gauss' theorem with the corresponding boundary conditions. To derive the formulation of Pa, we define a 
color function (i.e., a sharp characteristic function) as: 


<W x ) = 


x £ 

x £ f2 a . 


(27) 


and its smooth counterpart as 


^ (x)= Hi 


(ft('x.')W (x — x 7 , h r )dx'. 


(28) 


The gradient of <j> can then be found from Eq. 28 as 

V0(x) = fff^')V x W( X -x',h r )dx'. 


ClL)Ci a 

Using the definition of the surface delta function [43]: 

<5[n(x s ) • (x s — x)] = n(x) • V</)(x), x G 9 U Q a , for x s £ T, 

and noting that 

lim d) = d), 
h r -> o 

we can rewrite the surface delta function in terms of 0 as: 

<5[n(x s ) • (x s — x)] = n(x) • ^lim^ V</>(x), for x £ £7 U f2 a , x s 


£ 


r a . 


(29) 


(30) 


(31) 


(32) 
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The surface integral can then be rewritten as a volume integral through the surface delta function: 

ff ap(x ' s , t)dx' s = fff o:p(x, £)<5[n(x s ) • (x s — x)]dx, for x s £ T a , 

T a f2Uf2 a 

-ill a.p r (x, t)n(x) • V0(x)dx. 


(33) 


To uniquely define P;q(x, t ), we require it to vanish at a normal distance greater than h r from T a and require that 

^lim 0 JJJ Pn(x, £)dx = JJ ap(x' s ,t)dx.' s . (34) 

f2uf} a r a 

Comparing Eqs. 33 and 34 yields an expression for Pn (x, t) as: 


P n (x, t) = o:p r (x, £)n(x) • V</>(x), for x £ 17 U 17 a . 
Eq. 24 can then be rewritten by combining Eqs. 35 and 29 as: 


(35) 


dp r (x, t ) 
dt 


= V • (D(x)[Vp r (x, £) + /3p r (x, £)VV^(x)]) 


ap r (x, £) fff n ( x ) ‘ [0(x 7 ) V x tt;(x — x 7 , /i r )]dx 7 , for x £ 17 U !7 a . (36) 


Since p r is not uniquely defined on 17 a , we introduce a one-sided formulation by approximating Eq. 36 as: 


dp r (x, t ) 
dt 


= V • (D(x)[Vp r (x, £) + /3p r (x, £)VV^(x)]) 


• ap r (x, £) fff[ n (x) + n(x / )] • [</>(x 7 ) V x it;(x — x 7 , A r )]dx 7 , for x £ 17, (37) 


subject to the reflective Neumann boundary condition (Eq. 10). Note that </> is non-zero only in !7 a , where it is 
equal to 1 as defined in Eq. 27. Thus, the modified governing equation takes its final form as Eq. 9. 

Appendix B Calculation of Reaction Rate 

Similar to the derivation in Eq. 33, using the definition of the surface delta function and given pbuik = 1, the 
expression for the reaction rate can be rewritten as 


k on — fff M*) • J(x, £)][n(x) • V0(x)]dx 


(38) 


Substituting Eq. 29 into the above equation and using Eq. 27, a new expression of k on can be obtained: 

fc °" = /// [n(x) ' J(xt)] /// n(x) • V x iu(x — x / , fa r )<±x 7 dx. (39) 

f2ufi a fi a 

Similar to Eq. 37, the corresponding one-sided formulation is: 

k 0 n = fff [ n ( x ) ‘ J( x ’ £)] fff [ n ( x ) + n ( x/ )] ' V x iu(x — x 7 , h r )dx!dx.. (40) 

f2 a 

If the Robin BC (Eq. 5) is enforced, Eq. 40 can be reduced to 

‘■■Iff ap r (x, t ) Ilh x) + n(x 7 )] • V x tt;(x — x 7 , fa r )dx 7 cZx. (41) 
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Table 1 Comparison of Debye-Hiickel fits vs. ionic strength between experiment and simulations. 
RMSD of data to experimental k on is calculated over the range of ionic strengths between 0 and 0.67 
M. The unit of k® Q , k^ n and RMSD is lO^M - 1 min — 1 . And the error is the standard deviation of 
parameter fits using nonlinear least squares. 
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Figure 2 Comparison of SPH solutions to the analytical solution for the Smoluchowski equation 
subject to the Dirichlet BC on both r = R± and r = R 2 at steady-state with the relative L 2 
errors for different resolutions. Specifically, L 2 = 0.0666 for Ax = 8A(green square), L 2 = 0.0321 
for Ax = 4A(blue circle), and L 2 = 0.0153 for Ax = 2A(red triangle). 
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Figure 4 Panel A shows the discretized domain with R 2 = 128 A and the AchE molecule in the 
center with the reactive boundary shown in purple. Light blue indicates the outer boundary ( R 2 ), 
blue the solvent, green the protein and magenta the first (outermost) reactive boundary. Panels B 
and C show reactive boundaries 1 and 4, respectively in magenta spheres. 
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Figure 5 Contour of concentration distribution around mAchE (shown in dark gray) with the 
Robin BC (q = 8x 10 3 ) on reactive boundary 1 at steady state with a range of ionic strengths. 
Reactive boundary 1 is shown in purple. 



C = 0.0347ps t = 0.0693ps t = 0.191ps 


Figure 6 Time evolution of the concentration distribution around mAchE (shown in gray) with a 
Robin BC (a = 8 x 10 3 ) on reactive boundary 1 at 0.15 M ionic strength. Reactive boundary 1 is 
shown in purple. 
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Figure 7 (Left) k on as a function oft on reactive boundary 1 at different ionic strengths. Black 
square: 0.05M; red right-pointing triangle: 0.10M; blue asterisk: 0.15M; green circle: 0.20M; 
magenta diamond: 0.50M; cyan triangle: 0.67M. (Right) k on as a function oft on reactive 
boundaries 1-4, respectively, at 0.15 M ionic strength. Black square: reactive boundary 1; red 
circle: reactive boundary 2; blue diamond: reactive boundary 3; green triangle: reactive boundary 4. 
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Figure 9 Reaction rates of mAChE on reactive boundary 1 obtained from different methods. 
Black: from experimental data [4] (symbol) and fitted (line) to the Debye-Hiickel limiting law 
(Eq. 1); blue: from BD [9]; red: from FEM with Dirichlet BC [6]; green: from SPH with Dirichlet 
BC; magenta: from SPH with Robin BC using a = 8 X 10 3 . For standardization, both computed 
and experimental data are fitted to the Debye-Hiickel limiting law. 
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Figure 10 Illustration of the simulation domain and all boundaries: T 5 indicates the outer 
boundary, r m the molecular surface, and r a the reactive boundary 1; indicates the problem 
domain between and T a U r m . 







